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Abstract 


In a three part series of papers, a generalized finite element analysis 
scheme is developed to handle the steady and transient response of moving/ 
rolling nonlinear viscoelastic structure. This paper considers the develop- 
ment of the moving/rolling element strategy, including the effects of large 
deformation kinematics and viscoelasticity modelled by fractional integro- 
dif ferential operators. To improve the solution strategy, a special hier- 
archical constraint procedure is developed for the case of steady rolling/ 
translating as well as a transient scheme involving the use of a Grunwaldian 
representation of the fractional operator. In the second and third parts 
of the paper, 3-D extensions are developed along with transient contact 
strategies enabling the handling of impacts with obstructions. Overall the 
various developments are benchmarked via comprehensive 2- and 3-D simula- 
tions. These are correlated with experimental data to define modelling 
capabilities . 
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I. Introduction 


While extensive effort has been given to handling the response of 

f 1-3] 

structure with fixed loading sites , much less is available for prob- 
lems wherein either the load or contact regions move ^ . In this 

context, such structures as tires, rollers, compressors, bearings, pumps, 
turbines and manufacturing equipment (lathes, grinders, mills), all are 
typically excited by some form of traveling/moving load or boundary condi- 
tion. Such problems have been partially addressed in the recent works of 
several authors ^ . For instance, Padovan et al. ^ introduced 

the use of so called moving coordinate systems to handle traveling/rolling/ 

rotating structural systems. This includes both total and updated formula - 
f 91 

tions 1 1 involving primarily Hookean and Kelvin Voigt type formulations 
and 2-D benchmarks ^ 

In the context of the foregoing this paper series will extend previous 
work by considering the development of large deformation viscoelastic FE 
formulations for steady and transient traveling/rolling/rotating structure. 
Overall special emphasis will be given to: 


i) The development of large deformation viscoelastic moving finite 

[13] 

element (FE) formulations involving nonlinear elasticity 

and damping characteristics modelled by fractional integrodiffer- 

ential operators 


ii) The creation of a specialized hierarchically constrained Newton 
Raphson solution scheme for steady rolling problems; 

iii) The development of an implicit transient solver which incorporates 
a Grunwaldian type simulation of the fractional integro- 
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differential operator used to simulate the viscoelastic material 
behavior; 

iv) The development of 3-D translating/rolling cubic isoparametric 
type elements; 

v) The development of rolling/ translating type shell elements; 

vi) The creation of 3-D moving contact strategies; 

vii) Comprehensive 3-D simulation of a steadily rolling tire where 
analysis is correlated with experimentally derived tire test 
data; 

viii) The creation of transient rolling contact schemes capable of 
dealing with impacts with holes/bumps; 

ix) Integrating rolling/translating FE scheme with transient contact 
methodology; and, 

x) Comprehensively benchmark scheme with transient simulations of 
tire roll over events involving holes and bumps. 

Due to the number of topics covered, the overall paper is separated 
into three parts. The first covers items i) - iv). The second outlines 
iv) - vii). Lastly, the third covers items viii) - x). 

II Governing Field Equations 

To enable the modelling of rolling/moving structure undergoing large 
deformations, we shall employ the 2nd Piola Kirchhoff stress and Lagrangian 
strain tensor combination Secondly, to handle the moving/rolling 

motion, a change of reference base is used to define the requisite recti- 
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linear, Coriolis and centripetal accelerations. This yields the following 
equations of motion, namely 


( V 6 i k 


u i,k ))>J 


f. = 
1 


d 

o — r (u.) + 2pe. — (u ) + pe. e ft,u 
p , 2 v l K lmn m dt n mr mnk r k m 

dt 


( 2 . 1 ) 


where 6.., S.., u., f., Q , o and d()/dt respectively define the Kronecker 
ij ij i i m 

delta, the second Piola-Kirchhoff stress tensor, the displacement, body 
force, rotation vector, initial density and comoving time derivative. 

Since viscoelastic media are considered, 



S i j ( L 1 1 * * " * ’ dt 


( 2 . 2 ) 


such that L. . the Lagrangian strain tensor is defined by 

ij 


L. . = -z (u. . + u. . + u .u .) 
ij 2 i>i J » i i £ * J 


(2.3) 


The boundary conditions associated with (2.1) - (2.3) take the form 


i) For x.e3R ; 

i a 


S (6 . , + u . ,)n. = S. 
jk ik i,k 3 i 


(2.4) 


ii) For x.e3R ; 

l u 

u. -U. (2.5) 

i l 

where 

9R * 9R + 9R (2.6) 

u a 

such that 3R defines the surface of R. Since we are considering translating/ 
rolling structures, the loads defined on certain of the boundaries are of 
the moving type namely 


= S i (x - c(t),t) 


(2.7) 
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such that the vector x defines the initial coordinates, c(t) the trajectory 
undertaken by the load and t is time. Based on the functional dependence 
defined by (2.7), the various dependent field variables have the following 
form that is 



= s.jCn.t) 

B u i (n,t) 


L. . 

ij 


Lij(rj,t) 


( 2 . 8 ) 


where 

n = n(x,r,e,t) (2.9) 

[9] 

with n denoting the Galilean type position vector , and (r,0) the radial 
and circumferential locations relative to the center of rotation. 

Employing the functional representation defined by (2.8), the comoving 
derivative takes the form 


(;> - It ( "> + T <;> 


( 2 . 10 ) 




(u) + 2T (— (u)) + y(u) 


( 2 . 11 ) 



( 2 . 12 ) 


where the differential operators V and y are defined in the Appendix. For 
steadily translating/rolling structure, (2.10) and (2.11) reduce to 


3t V - , ' ( S ) 

and 

A 1 

^“7 (u) » y(u) 
dt ' 


(2.13) 


(2.14) 


To enable the derivation of the requisite FE formulation, (2.1) - 
(2.6), (2.10) and (2.11) can be employed to yield the virtual work prin- 
ciple cast in moving coordinates. After several manipulations we yield 
the expression 

2 

/ {6(L) T S + 6(u) T ( [ m. ] (u) + 

R ~ 3t Z ‘ 

[[m^] + 2[m 1 ]'{'( )] (u) + [[m^] + 

[m 2 ]'F( ) + [m 1 3 Y ( )]u)} dv 


f S., (6,. + u, . )n.du,ds 


3R 


jk v "ik ■ “i,k /w j““r 


(2.15) 


where ( ) denotes matrix transposition, 6( ) the first variation and 


T 

S = 

^ S ir S 22 ,S 33 ,S 12 ,S 23 ,S 13^ 

(2.16) 

II 

H 
t-3 i 

^ L 11 ,L 22 ,L 33 ,L 12 ,L 23 ,L 13^ 

(2.17) 

[nij] 

5 p[I] 

(2.18) 

[m 2 J 

= 2pe. n 
lmn m 

(2.19) 

[m 3 ] 

= pe , e Q, 

inr mnk r k 

(2.20) 


For the case of steady rolling/translating, (2.14) reduces to the form 


/ {<5(L) T S + 6(u;[ [m- ] + [nuM ) + [m.M )]u} dv = 
R " 

/ S 41 .(6.,. + u. u ) n^du^.ds 

3R 


jk lk i,k' 3 i 


( 2 . 21 ) 


[9] 

where here due to the use of Galilean type coordinates » time is removed 


from the formulation for both linear and nonlinear situations. 
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The preceeding formulations were cast in so called total Lagrangian 

form wherein all kinematic fields are referenced relative to the initial 

configuration. Alternatively, for very large deformation situations, the 

updated Lagrangian strategy can be used. For such a formulation, the 

reference base is reassigned intermittantly to intermediate kinematic 

states. In this context u., S.. and L. . are recast in the form 

i ij ij 


u. = V. + u. 

1 1X1 

L. . ■ L. . + L. . (2.22) 

ij ij t ij 

S. . ■ S. . + S. . 

IJ IJ T IJ 

where (V - ., L.., S..) and (u., L.,, S..) are respectively the dis- 

1 ij’ ij X 1 x ij x ij r 1 

placement, strain and stress fields defining the reference state at x and 

their associated excursions. Based on (2.21), (2.14) takes the form 

2 

/ { 6 ( L) T S + 6( L*) T S + 6( u)([m.] ( u ) + 

R x~ x - T " " T ~ 1 3t 2 T ~ 


[[m 2 ] + 2[m 1 ]t[l] ( T u) + [[m 3 ) + 


[m 2 J'F( ) + [m 1 ]Y( )J t u)} dv = 


/ ( S„ (6, u + V, + u. J + 


a r 

x 


X jk v ik i,k x i , k / 


S .. u. . ) n.6( u. ) ds 
jk x i,k' j t i 


where 


6( t L) 5 ± {«( T » lfj + 4^^) + 


v t ,i a( T u i .j ) + V M S< t u ll ,i ) + 


t u e,i 5( T u t,i ) + i( i u i,i ) T u t.j ) 


(2.23) 


(2.24) 



For the steady case, (2.23) reduces to the form 


/ { 6 ( L) T S + 6( L*) T S + 6( u ) [[m_] + 

T~ T- T~ ~ T~ J 

K 
T 

[m 2 m ) + [m 1 ]Y< )] T u)) dv = 


/(_S..(6.. + V.. . + u, J + 


3 R 

T 


x jk v ik i,k x i,k 


S .. u. ,)n.6(u.)ds 
jk x i,k j v x i' 


(2.26) 


III. Viscoelastic Material Properties 


To complete the set of governing field equations, constitutive rela- 
tions need to be defined. For the current purposes we shall consider the 
translating/rolling of viscoelastic structure. Since rolling/translating 
on smooth or uneven surfaces involves a wide range of speeds and potential 
spectral excitations, the viscoelastic formulation must be able to accommo- 
date such factors. In this context it is well known that the use of the 
traditional Kelvin Voigt simulation cannot handle problems involving a 

perfusion of spectra. Because of this we turn to the use of fractional 
integrodif ferential representations . As shown by Bagley and Torvik, 

such operators have an improved capacity to handle frequency varia- 


tions. 

To enable the development of such a formulation, we recall the classical 
generalized Maxwell-Kelvin-Voigt type representation, namely 


S ij + J 4j l ab 71 < S ab> ' G ij (L ir L 22--> + 1 ^ <L -» ) 

uU 


i u* dt* ab 


(a) ( o) 

where he^e ” . and a , are respectively the compliance and 

■"ijab ijaD 
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stiffness properties of the various stress and strain rate terms. Note 

tvior 
[13] 


G. .(L, , , . . . ) defines the nonlinear elastic behavior which say in the case of 
ij 11 


rubber would be the Mooney-Rivilin correlation 


Reinterpreting the 


integer ordered differential operators appearing in (3.1) to be of the more 
comprehensive fraction type ^ , we yield the expression 

S. . + E u^.D (S , ) = 
ij l ijab ab 


G ij (L ll ,L 22 ,,,,L 13 ) + £ <x ijib D q £ (L ab ) +P ijab dt (L ab ) 
such that here in Riemann-Liouville form 


(3.2) 


V s *’ ■ C Zm 


P„ ab 


with T defining the gamma function 


S ab ( T) dx 
[18] 


(3.3) 


Alternatively, D ( ) can be 


expressed in the computationally more convenient but equivalent form 


developed by Grunwald 


[16] 


namely 


D (S) = E (At)" P * A(p 0 , j+l)S(t-j At) 


(3.4) 


such that 




A(p £ ,j+i) (-D J (j ) - r(-p £ )r(j+i) 


(3.5) 


Note, the coefficients A(p £ ,j+1) denote the memory of the material. For 
appropriately chosen p £ , they represent an essentially monotone decreasing 
set. As p £ the order of the fractional derivative can be experimentally 
chosen [1^*15]^ a w ^<j e ranging series of histories can be accommodated by 
(3.2). 

To expedite the development of the requisite solution algorithm, (3.2) 


is converted to matrix form namely 


9 


s + E [u (2,) ]D (S) = 

l P n 

G(L) + 2 [a U) ]D q (L) + [g] ^ (L) 

& 2 . 


(3.6) 


where here 


r £, _ U) 

1)1 1 = w ijab 

,JO, = ft) 

1 - “ijab 
[S] S e ijab 


and based on (3.4) we see that 


D (S) = (At)' P d E A(p , ,.)S(x,t-jAt) 
P Z ~ j=0 *0+1 - - 


(3.7) 


(3.8) 


D (L) = (At) E A(q 1 )L(x,t-jAt) (3.9) 

q SL ' j=0 *’ J 1 " " 

Since we are dealing with rolling/translating structure, recalling the use 
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of the Galilean transform it follows that (3.6) converts to the form 


S(n,t) + 

E [y (ll) ](At)" P * E A(p £ ^ + 1 )S(n(x-c(t-jAt)),t-jAt) = 


G(L(n,t)) + [p](Y(L(n,t)) + ^r (L(n,t)) + 


E 

i 


[a^](At)" q £ E A(q £ ^. + 1 )L(n(x-c(t-jAt)),t-jAt) 


(3.10) 


To recast (3.10) into a more convenient form, it must be recognized 
that the family of terms defined by S(n(x-c(t-jAt)) , t-jAt) and 
L(n(x-c(t-jAt)) ,t-jAt) je[0,°°) represents the history of a given particle. 
Obviously the particle occupies different locations during its history. For 
simplicity we shall designate the succession of positions by the nomenclature 
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(r,m,t+At) such that r designates the radial position, m the circumferential 
positions and t+At the time. 

For rotating bodies, the m t ^ 1 circumferential location is defined by 
jAn r wherein 

An * c it (3.11) 

-r -r 

In this context at time t the location of the particle is given by 

(r m) 

(r,m-l,t). Based on such nomenclature, S is cast in the form t+At S * 

S (r ’ m > = s(n + (m+kJ )An , t+At) (3.12) 

t+At- - ~r r ~r 

such that n r is the initial position, k the total number of complete revolu- 
tions, the number of time increments required to define a given revolution 
and Ari r the incremental circumferential motion of the r ^ position. In a 
similar context, we see that 


t+At L (r ’ m) = L(n„ + (m+kJ^)An^,t+At) 


- -r 


r -r 


(3.13) 


Employing the foregoing nomenclature, it follows that (3.6) reduces to 


t «t5 (r '" ) + l lvWm ei ( t+it5 <r ’ m)) - 


- ( r+At- (r ’ m)) + [P]( ' F( t+At- (r,m)) + 


d t r(r,m)\ . .. r«^) in ( j(r,®)\\ 

at ( t+At- ) l l 1 q/t+At^ 


(3.14) 


Generally the solution to FE simulations involves the use of tangent 
type material formulations. In this context, we note that 


t+At- 


(r,m) = Ag (r,m) + g(r,m-l) 


t- 


(3.15) 
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t+At- 


(r,m) _ AL (r,m) + L (r,m-1) 


and 


G( L (r>m) ) - [Di r -" , J AL (r > m) + G(.L (r>m ' 1) ) 
~ t+At T - - t~ 


(3.16) 


(3.17) 


where here [D^] defines the tangent elastic properties. Based on these 
expressions, it follows that (3.14) reduces to the form 


AS 


(r,m) + g(r,m-l) + 


t- 


Z [y (£) ](At)" P a {A(p ,l)(AS (r,m) + S (r,m ' 1) ) + 

£ 16 

I A(pl.j+l) t . (J . 1)at S (r - n - j) ) - 

J=1 J 

G( t L (r » m ' 1) ) + [D^ r,m) ] AL (r,m) + 

Z [a (£ *](At)‘ q £ {A( qjl ,l)(AL (r,m) + t L (r » m_1) ) + 

1=1 

[p] {'F(AL (r,m) ) + 'P(.L (r,m " 1) ) + 


|- (AL (r - m) + ,L (r ’ n, - 1) )} 

ot tl ~ 


(3.18) 


At time t the prior location of the particle is (r,m-l), hence we have that 

g(r,m-l) + 


Z[y (0 ](At)" P * Z A(p ,j+l) S 
£ j=0 J 


(r,m-l-j) _ 


G( t L (r ’ m ‘ 1) ) + 


(over) 
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Z [a U) ](At)' P * Z A(q ,j+l) L (r » m " 1_j) + 

% j=0 1 1 JAt ~ 

[p] {'F( t L (r,m_1) ) + §£ ( t L (r ‘ m_1) )} 

Differencing (3.18) and (3.19) we yield the incremental expression 

[y (r ’ m) ] AS (r ,m) = 

[a^ r,m) ] AL r,m) + [p] 'F(AL ( ' r,m) ) + 

r o i £_ f A r (r » m) i + h ( r » m ) 
m 3t (AL ) + t+At h ya(3 

where 

[y (r,m) ] = [[I] + Z (At)' p * A(p ,i)]y (4) ]] 
l 

[a (r,m) ] = [[D^ r,m) ] + Z (At)" q H A(q.,l)ta (£) ]] 

1 l * 

and 

t+4t~^yo™ ) * \ (lt) ' qi! 1“ U>1 } l 0 - <1- 

sgn(J,0))A(q t , j+l)) t _j 4t L* r,m 1 ^ - 

Z (At)" p a [y (l) ] Z (A(p.,j+2) - (1 - 
l j=0 

(r,m-l-j ) 


(3.19) 


(3.20) 

(3.21) 

(3.22) 


sgn(j,0))A(p 4 ,j+l)) t _j At S 


(3.23) 
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IV. FE Formulation 

For the current purposes, a displacement type representation is employed 

[ 2 ] 

to develop the requisite FE formulation. In this context we have that 

u = [N]Y (4.1) 


where [N] is the shape function and Y the nodal deflections. For transient 
situations it follows that 


and 


[N] = [N(n>] 

[ 10 ] 


Y = Y(t) 


(4.2) 


(4.3) 


Based on (4.1), (2.10) and (2.11) yield the following expressions for 
velocity and acceleration, namely 

^ (u) - V ( CNJ )Y + [N] (Y) (4.4) 

2 2 
^ (u) = y ( [N] )Y + 2Y([N]) fr (Y) + [N] ^ (Y) (4.5) 

dt 2 “ " at • at 2 ~ 

[9] 

Continuing, 6L takes the form 

6L = [ B* ] 6 Y (4.6) 

Employing (4.4) - (4.6), (2.14) can be used to develop the following FE 
formulation, that is 

/ UB*] T S + [N] T ( [m. ] [N]Y + 

R " 

[[m 2 ][N] + 2[m 1 ]'P([N])]Y + [[m 3 ][N] + [n^MfN]) + 

[m 1 ] Y ([N])]Y)} dv = F (4.7) 
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In the case of steady translating/rolling Eq. (4.7) reduces to the form 

/ {[B*]S + [N] T ([m~][N] + 

R 

[m 2 M[N]) + } dv = F (4.8) 

Since (4.7) and (4.8) are inherently nonlinear, numerical 
schemes are required to generate their solution. This will be developed 
in the next section. 

V. Solution Algorithms 

The FE simulations given by (4.7) and (4.8) respectively define transient 
and steady state aspects of the translating/rolling problem. Solution algo- 
rithms to such models are developed in the succeeding subsections. 

V. 1 Steady Problem 

Since (4.8) is inherently nonlinear, the Newton-Raphson scheme will be 
employed to affect its solution. To improve the resulting algorithms* 
stability and self adaptiveness, a locally constrained version will be devel- 
oped. To start, (4.8) is recast in incremental form namely 

/ {[G] T [S][G] + [B*]AS + [[m 3 HN] + [m 2 M(N]) + 

[m 1 ]y([N])]&Y } dv = AF (5.1) 

such that [S] is the prestress matrix. At this stage, it must be recognized 
that due to the viscoelastic properties, (5.1) must reflect the appropriate 
history effects on AS. In this context, for each point in the domain given 
by R, AS and [S] must be linked to its past so as to define the requisite 
history effects. Hence, (5.1) is recast as follows: 
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/ {[G (r ’ m) ] T [ t+At S (r>m) ][G (r ' m) ] + 


[B* (r ,m) ] AS (r,m) + [[m,][N] + [m,]?([N]) + 


[m 1 ]Y([N])]AY (r,m) } = AF (r * m) 


Employing (3.20) in conjunction with (5.1) we yield the following 
assembled algorithmic expression 


. AY = . , [K ] ~ 1 ( .AF - .H Q ) 
l - i-l l l - i~yap 


i-lIKl = i.jUKj) + [K pp ] + [K p ]] 


where 


i-itV = 1 + i _ 1 ([B*] T [y]” 1 [a][B*])} dv 


i-l [K yp ] = 1 [B*] )) dv 

R 


. ,[K ] = /[N] T . , ( [m_] [N] + [m-M[N]) + [».h([N])) dv 

P - 1" 1 J Z 1 


(5.2) 

(5.3) 

(5.4) 

(5.5) 

(5.6) 


.AF = F - / . . { [B*] T S + [N] T ( [m_] [N] + 
1 " R 1-1 


and 


[m 2 ]'F([N]) + [m 1 ] Y ([N]))Y} dv 


,-1 


(5.7) 


(5.8) 


.H a - S [y] .h 0 dv 
i~yap ^ 1 -yap 

(r m) 

such that [y], (a] and h^ are the assembled versions of [y v * '] , 

[a (r - n) ] and h (r> " ) . 

- yap 

Note, for rolling geometries which are symmetric, the various mechanical 
fields are periodic. Because of this, the computation of „ is greatly 
simplified. Specifically, turning to the instantaneous local (r,m,t+At) 
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particle location, we see that due to periodicity, for the previous roll 
the (r,m-J r ,t+At-J r At) position possesses the same mechanical fields. 

In this context, it follows that 

f . A J (r, ” ) = S(n + ( m+k J ) Ah , t+A t ) 

t+At~ ~ ~r r ~r 

= S(n r + (m-l+kJ r )Ar) r »t+(l-J r )At) 

Hence, for the steady rolling case we have the following interrelationship 
between succeeding rolling states, that is 

t+At? (r ’ m) = ? (q r + (">+kJ r )An r ,t+At) 

s (r,m) = + (^j ) An ) 

~ ~r r -r 

s (r,m) = S ( n ^ + mAr ^) ( 5 . 9 ) 

for all ras[l,J ] and ke[0,°°). Similarly, 


L(n r + (nH-kJ r )Ari r ) = L(n r + mAr| r ) (5.10) 

Based on the mechanical history defined by S^ r,m ^and (3.23) reduces 


to the more convenient form 


h (r, J } h Z (At)' q «- [a U) ] Z C(q ,j,m,r) L (r ’ j) - 
~ yap l j-1 £ 

J r 

Z (At)" P * [y ( °] Z C(p , j ,m,r) S (r,j) 


(5.11) 


where 


C(q £ , j »m, J f ) * 2 (A(q^,m+l-j + (sgn (j,m) + k)J r ) - (1- 


sgn ( j ,m) )A(q^,m-j + sgn (j,m) + k)J r ) 


(5.12) 
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C(p , j ,m, J„) = 2 (A(p^,m+1- j + (sgn (j,m) + k) J r > - (1 - 
* ‘ k 

sgn (j ,m))A)p Ji ,m-j + (sgn (j,m) + k) J r >) (5.13) 


The coefficient families C(q £ , ...), CKp^,...) represent the history 
effects on the (r,m) th particle. For instance, (5.12) defines the influence 
on the strain state of the succession of kef l,® 0 ) rolls through the (r,j) 
orientation. Similarly, (5.13) defines the influence on the stress state 
of the succession of ke[l,°°) rolls through the same orientation. As seen 
from (5.11), the net (r,m) t * 1 history is obtained by summing through all the 
je[l,J r ] positional effects for all ke[0,«>) rolls. 

The algorithm defined by (5.3) yields the succession of ^AY converging 
to the solution 


Y. 

~i 



+ .AY 

i ~ 


i=l,2,3. . . 


(5.14) 


To improve the robustness of the overall solution process, local constraints 
are introduced. Since rolling structures have a variety of substructural 
zones undergoing different levels of nonlinearity, AY is controlled separate 
ly through the use of a partitioned constraint process. Namely, (5.3) is 
replaced by 

.AY = diag [X . *(F - .H yajJ ) - .^[kf 1 f [B*] T S)dv (5.15) 

R 

where diag [A^] is defined such that ^AY is hierarchically controlled. In 
particular 

diag [A i J ■ 


V 1 ! 1 


A i2 [I 2 ] 0 


X. [I ] 

IK K 


0 


(5.16) 
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such that A.„, £e[l, ] are the various constraints and [ 1 0 ] ; £.e[1,k] a 
set of identity matrices whose individual sizes define the groups of 
degrees of freedom controlled by each given A^. The choice of the A_^; 
£e[1,k] family is obtained by either: 


i) Requiring each of the various partitions of ^AY to remain upper 
bounded by a specified linear constraint; or 


ii) 


By employing a more general functional constraint say as in the 

[19,20] 

hyperelliptic constraint surface of Padovan et al. * 


Before establishing either of the foregoing schemes, all the various 
key vectors and matrices must be appropriately partitioned. In particular 
we let 





A T 

i^2* 


.Ay T ) 

1 J -K 


(5.17) 


( i-l [K] " 1)T = t ± tE l 1 * i tE 2 ] 


i^c” 


(5.18) 


with 


i V E 




‘W + 


'V 11 


(5.19) 


such that the various partitions ^y^ and ^[E^], £e[1,k] have respectively 
the same number of rows as the identity matrices [I^]; £e[ 1 »<]» For sim- 
plicity, the driving potential F - is recast as 


.F = F - 


•H a 

l yap 


(5.20) 


Hence, employing (5.17), (5.18), and (5.20), (5.3) reduces to 


*,T. 


± Ay a * - i [E £ ]/ i-l ([B ] '5 )dv ; tel 1 ’*! 


(5.21) 


[19] 

Considering the elliptic constraint function approach , we intro- 


duce the localized version 
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where £e[l,<] and 

iU m i-A * i iy -i 


(5.22) 


(5.23) 


ill - (A i£>i? 


(5.24) 


such that defines the aspect ratio of the ellipse. Its choice will be 
discussed later. Solving (5.21) - (5.24) simultaneously, we yield a quad- 
rative expression in namely 


i A £l^ A i£^ + i A £2^ A i£^ + i A £3 


(5.25) 


Since the coefficients are known scalar s, A^ can be directly 

evaluated from (5.25). Note, the roots chosen will be dependent on the 


state of definiteness of [k_ q ]. 

TyBp 

The choice of the size of the [1^] partitions and their associated 
warping parameters 0^ can either be user defined or automatically/self - 
adaptively generated. If user chosen, the size of the various partitions 
can be taken from a variety of factors namely 


i) Inherent substructuring of a given system into component parts 
defined by 

• geometry 

• material groups 

• boundary conditions 

or 

ii) By geometry, material or boundary condition induced nonlinearity. 

Such an approach also points to the appropriate choice of the local 
warps Note, the control of the size of allowable local field excursions 
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is obtained by requiring that they remain bounded by the family of elliptic 
constraints defined by (5.22 - 5.24). 

For instance, from blueprint information clearances are often given 
for various substructural components. Such information can be used to 
determine upper bounds on a given AY^ allowable. In terms of (5.22), such 
allowables would yield the following choice of 9^ namely 

•»i Mo^M^lhtallo-ablJI 2 (5 ' 26) 

To self adaptively update such partitioning and warping, running checks 
can be made on: 

i) The normed partition level out of balance loads; 

ii) The level of strain energy stored/dissipated; or, 

iii) The normed excursions in the associated mechanical fields. 

Once various combinations of i) - iii) are triggered, the appropriate parti- 
tions may be enlarged or shrunken along with associated changes in 0^. 

Noting the architecture of the constrained algorithm defined by (5.21) - 
(5.25), such partition and warping updates involve only straightforward 
accounting changes to define/locate a given partition. 

V.2. Transient Problems 

To handle transient problems, we shall adopt the implicit manner of 
formulation. In this context we recast (4.7) in so called incremental 
format, that is 

/ {[G] T [S][G] + [B*] T AS + [N] T ( [m. ] [N]AY + [[m ](N] + 

R ~ 

2[ mi M[N])]AY + [[m 3 ][N] + [m 2 M[N]) + [m^yC [N] ) ]AY)} dv = AF 


(5.27) 
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For the current purposes, the solution to (5.27) is obtained by employing 
the Newmark Beta method as defined by the following time stepping algorithm 
namely 


t+At? " a 0 A - a 2 t- " a 3 t- 


t+At- " J + a 0 t? + a 7 t+At? 


(5.28) 

(5.29) 


a o = 

l/(3At 2 ) 

a 4 = 

y/B-1 

a l = 

Y/(B At) 

a 5 

At/2(y/ 3“2) 

a 2 = 

1/ ( pAt) 

a 6 * 

At(l-y) 

a 3 = 

1/(23) - 1 

a 7 = 

yAt 


(5.30) 


Rearranging (5.28) and (5.29) into incremental form we obtain the 
expressions 

(5.31) 

(5.32) 


A Y = V? +C 1 J + C 2 t? 


AY = C-AY + C. Y + C c . Y 
3 - 4 g~ 5 t- 


where 


C 0 a 7 a 0 

C 1 = ” a 7 a 2 

C 2 = a 0’ a 7^ a l _1 ^ 


C 3 = a 0 

C 4 = ~ a 2 
C 5 = -(a 3 -l) 


Based on (5.31), (3.20) takes on the following form 
Ag (r,m) _ ^(r,m)j-lj|. Qt (r,m)jj B x (r,m)j + 

[3W[B* (r,m) ]) + C 0 [8][B* (r,m) ]]AY (r,,n) + 

r (r.m),-! (r,m) 

^ ] t+At- yap 


(5.33) 


(5.34) 


such that 
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(r,m) m 
t+At2 yotg 


t«ti < S^ ) + teHB' ,(r, ” ) ) {Cj + c 2 t Y <r '"- 1) ) (5.35) 

At this stage, employing (5.31), (5.32) and (5.35), it follows that at 
the fully assembled level (5.27) yields the following expression 


t+4t tIK T I + ‘V + [K P 114Y ’ t +1 t ( ?pae + 4 ? + Ip* 


where here 


r _ , - fii (r,m) 1 (r,m) 

t+At G yag £ t+At? yag 

K r 


AF = F - /[B X 1 T S dv I 
t+At - t+At- J - * 


t+At [K yg ] = f f B *l T fV‘]' 1 [C 0 [B][B 5!t ] + [g]T([B*])] dv 
R 


t+At 1 p 


[K ] = / [N] T [C 3 [m 1 ][N] + C 0 ([m 2 J[N] + 


R 


2[m 2 J'F([N])) + [m 3 ](N] + [n^lVUN]) + [mjMtN])] dv 


and lastly 

t+At I * / [N] T {[C 4 Cm 1 lEN] + C 1 ([m 2 ][Nj + 
R 

2[m 1 M[N]))] t Y + [C 5 [ mi ][N] + 

C 2 ([m 2 J[N] + 2[m 1 ]'P([N]))] t Y} dv 


Solving (5.36) for Y yields 

iY * t +1 t [K D I ‘ 1 t +4 t<5ya 8 + + V 


(5.36) 

(5.37) 

(5.38) 

(5.39) 


(5. AO) 


(5.41) 


(5.42) 


where the dynamic stiffness is given by the relation 
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t + 4 t [l V ‘ t +i t t[K T ] + (1 W + 1 k p ]1 

Based on (5.41), the time marching scheme defined yields 
solution 


(5.43) 
the overall 


Y 

t+Ati 


Y + AY 


(5.44) 


To improve the solution, (5.42) can be iterated in much the same manner as 
the steady case. Convergence for such a process is defined by the usual 
displacement norm check . 


VI. Discussion 

In the preceeding sections, a generalized nonlinear viscoelastic FE 
formulation was derived for translating and rotating structure. This in- 
cluded the development of both steady and transient solution algorithms. In 
the case of the steady formulation, a new locally constrained solution algo- 
rithm was created. For the transient case, an implicit algorithm involving 
the use of the Newmark Beta method ^ was developed. Note to enhance the 
capabilities of the overall formulation, the more comprehensive fractional 
operator is used to represent the viscoelastic effects. 

In parts II and III of this paper, the foregoing algorithms will be 
extensively benchmarked for both steady and transient situations. As noted 
earlier, this will include evaluating steady 3-D rolling contact simulations 
of pressurized tires. For the transient case rolling/translating over bumps 
and holes will be simulated. To enable such simulations, part II will 
include: 

i) The development of 3-D rolling/translating elements handling 
nonlinear kinematics and material properties; 
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ii) The development of rolling/translating shell elements; 

iii) The generalization of 3-D contact strategies to handle rolling 

contact; and, 

iv) Comprehensive 3-D benchmarks with experimental verification. 

In part III, the main thrust will be to: 

i) Develop transient rolling contact schemes capable of handling 
impacts with obstructions; 

ii) Develop the overall solution architecture integrating both the 

rolling/translating FE scheme with the transient contact scheme; 
and, 

iii) To comprehensively benchmark the scheme. 

Note the transient and steady translating/rolling FE simulations and 
associated solution algorithms have architectures essentially compatible 
with traditional formulations. Apart from the nonsymmetry induced in the 
static and dynamic stiffness, the changes required to incorporate such schemes 
in GP codes are minimal. 
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Appendix 

In section II the governing field equations were developed. Here the 
so called moving total/updated Lagrangian coordinates were introduced. 
Assuming a cylindrical type coordinate system to define the geometry of 
rolling/ translating structures, n can be defined by the expression 


n 


r cos (9+<(»(t)) 
r sin ( 0+d> ( t ) ) 


(A. 1) 


where here for constant rolling speeds about the x^ axis 
(J> (t) = ft t 


(A. 2) 


Based on (A.l) and the fact that u^ = 
total derivative is given by 


j . an. 

d? <u i ) * 5T + ST 'V at 

J 


u^(n>t), we see that the first 


(A. 3) 


where 


an, 

= - flr sin (0+<J>) 


(A. A) 


8n ? 

* nr cos (0+<|») 


(A. 5) 


Combining (A. 3) - (A. 5) we see that 


3t (u i> ' fe (u i> + * (u i ) 


(A. 6) 


such that the ¥( ) operator is defined by the expression 


¥( ) = -nr sin (9+<J>) f — () + nr cos (0+<|>) () 

dn^ or i2 


(A. 7) 
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Continuing, the second total derivative takes the form 


(u.) * 2-r (u.) + 2'F (|- (u.» + y(u . ) 

dt 2 1 at 2 1 at 1 1 


(A. 8) 


where here 


y( ) ■ -n 2 r cos (0+<t) |^j— ()-n 2 r sin (0+<J>) |^~ ( ) + 
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ani 


n 2 r 2 sin 2 ( ©+d» ) ^—r ( ) + fi 2 r“ cos 2 (0+$) T ( ) 

an 2 

2 2 3 2 

2ft r sin (0+<|>) cos (0+<J>) a _ a '~ ( ) 

& n ^ ® ’>2 


(A. 9) 


For the special case of steady rolling, it follows that - u^(q). 
In this context (A. 6) and (A. 8) reduce to the expressions 


dt <“i> * «»i> 


— - (u i ) - yCuJ 
dt 


(A. 10) 
(A. 11) 


where here time appears implicitly through the definition of n* As can be 
seen through (A. 10) and (A. 11) the time derivative is replaced by a spatial 


one. 
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Abstract 


In a three-part series of papers, a generalized finite element 
solution strategy is developed to handle traveling load problems in 
rolling, moving and rotating structure. The main thrust of this section 
consists of the development of 3-D and shell type moving elements. In 
conjunction with this work, a compatible 3-D contact strategy is also 
developed. Based on these modelling capabilities, extensive analytical 
and experimental benchmarking is presented. Such testing includes 
traveling loads in rotating structure as well as, low and high speed 
rolling contact involving standing wave type response behavior. These 
point to the excellent modelling capabilities of moving element strategies. 


i 


I. Introduction 


This is the second paper in a series of three considering the 
development of large deformation viscoelastic FE formulations for 
steady and transient travelling/rolling/rotating structure. Overall, 
this part will give special emphasis to: 

i) The development of 3-D translating/rolling 
isoparametric type elements; 

ii) The development of rolling/translating type shell 
elements; 

iii) The creation of 3-D moving contact strategies, and 

iv) The comprehensive 3-D simulation of a steadily rolling 
tire; here the analysis is correlated with experimentally 
derived tire test data to provide real world corroboration. 

Since part one of this series has provided a fairly thorough review 
of previous work, for the sake of conciseness, we shall immediately get 
into the development. Note, the overall structure of the paper follows 
that defined by items i) - iv) noted earlier. Of particular importance 
is the benchmarking phase which involves both analytical comparisons, 
as well as, experimental test data. The empirical numerical correlations 
include: 

i) Standing contact; 

ii) Frequency properties as per small deformation superposed 
on large and; 

iii) Full rolling contact through all possible ranges of 
velocity. 
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2. Three Dimensional Formulation 

Recalling Part 1 of this series, FE equations were derived for 
moving load problems involving: 

i) Transient inertial effects; 

ii) Large deformation kinematics; and 

iii) Viscoelastic properties. 

Based on the use of a transient version of the Galilean transform [1-3], 
the following moving formulation was developed. 

/ {[B*] T S + [N] T ([m ][N] Y + [ [nu] [N] + 

R 

2[n» 1 ] Y([N])] Y + [[m 3 ][N] + [m 2 ] Y([N]) + 

[m 1 J Y([N])]Y)}dv = F (2.1) 

For the case of steady state mot ion , the use of the multiply 
constrained partitioned Newton Raphson scheme [3,4] yields the fol- 
lowing solution algorithm namely 

J\Y = diag [A ± ] ^[K] 1 (?- i H j ja p) ' 

i-i ^ 1 " 1 1 i-i (tB * ]T ? )dv (2,2) 

R 

For this case, the consistent mass matrices are embedded within [K], 
that is 

i-i [K1 - 1-1"%' + [ V + (k p" (2 - 3) 


where 
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i -1 


[ V 


/ [N] i _ 1 ( [m 3 ] [N] + [m 2 W[N]) + [n » 1 J't ( [N] ))dv 


(2.4) 


To establish . . [K ] for 3-D formulations, we shall employ a 20 node 
l-l p 

isoparametric serendipity type brick element. Noting Fig. 2.1, it has 
three displacement degrees of freedom per node. The displacement fields 
in the 20 node solid element depend quadratically on the position within 
the element [5,6]. That is, the components of the shape function of this 
element have quadratic terms involving the isoparametric coordinates. Fi- 
gure 2.2 shows the element in isoparametric space. The components of the 
shape function of this element are given in Table 2.1. 

Noting the differential operators ¥( ) and y( ) given in (2.4), the 
shape functions must be differentiated spatially. Because isoparametric 
elements are used, the first step is to relate the derivatives with res- 
pect to the global Cartesian coordinates (x^, x^, x^) to derivatives with 
respect to the elements isoparametric coordinates (€^» £ 3 )* This is 

necessary because the shape functions of isoparametric elements are 
written with respect to the £ coordinates. By convention, these shape 
functions are used to relate the displacement fields in the element to the 
nodal displacements. 

Noting (2.4), it is seen that expressions for 


{ 



3(n^, n 2 » n^) 





t 


(2.5) 


4 


and 


r 3 2 u. 

i 


3n, 


3 2 u. 

1 


9(n x , n 2 , n 3 )^ 


- 


3n_ 2 


3 2 u. 

i 


3n, 3n_ 
1 J 


( 2 . 6 ) 


must be found in terms of the nodal displacements. The forms used in 
(2.5) and (2.6) are shorthand for the full number of first and second 
derivatives, namely 

(3u^/3n. , 3 2 u./3n i 3n^); (i,j)e[l,3]. 


The first derivatives are related by 


3u. 3u. 

Ufr- h — r-T> = [J] (— — - — — } 


Hty C 2 , C 3 ) J 


3(n^, n 2 » n^) 


(2.7) 


where [J] is the Jacobian matrix given by 


[J] 


t — — ] I 1 ” r ° W 
l 3Cj J l j ■ column 


( 2 . 8 ) 


This expression can be inverted to give the relationship between the 
Cartesian and isoparametric spaces needed for the first derivatives, 
that is 


3u j 

^3(n^, n 2 » n^)^ ^ ^ ^3(5^, C 2 > £ 3 ) 


} 


(2.9) 
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To establish the requisite second derivatives, the chain rule is 

2 2 

applied twice. In the context of the 3( )/5 3 and 3 ( )/d^ 1 derivatives, 
such an operation yields 
3u. 


3 3u. 3n. 

= *=1 ^ ^ 


( 2 . 10 ) 


and 


3 2 u, 


se j 


= Z 




“ 35. '3n„ 35, 


d=l 


( 2 . 11 ) 


which upon expansion and combination of like terms reduces to 

a2u i 1 a v ^ *\ , .. a2 "i + 

IT 5 ' t-i ‘an 2 S S "* *J 3 "i 3 "2 S «j 35 J 

J ^ 


2 2 
3 u. 3n. 3n 0 3 u. 3n 0 3n- 

2 i i — 2. + ? 1 

3n 1 3n ;J 35.. 3S.. 3n 2 3n 3 35.. 3^ 


( 2 . 12 ) 


Repeating this process for the remaining second derivatives, we obtain 
the following overall operator expression namely 


3 2 u. 


3 2 u. 


{ 3(5 1 ,C 2 .e 3 ) 2} * tJlJ {a( V n 2’ n 3 )} + U21 { 3(n 1 ,n 2 ,n 3 ) } 

The matrix [J^] is (6x6) and [J 2 1 is (6x3). Now, based on (2.9), (2.13) 


can be inverted to yield 


a 2 u. 


3(n 3 »n 2 »n 3 ) 


,} = [J 


, 3 2 u. 3u. 

J <— “ — ; 2 >* [J 3 ] { 3u^l } (2 - 14) 

a(5 r 5 2 ,c 3 ) 12 3 


where 
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[J 3 1 » [Jj 1 t J 2 ] 


(2.15) 


The first and second derivatives on the right side of (2.14) are 
with respect to the element's isoparametric coordinates. The shape 
function and hence the displacement fields are described using this 
coordinate system. Recalling Part 1, the displacement field in the 
element is related to the nodal degrees of freedom by 

u = [N]Y (2.16) 


For the components, (2.16) yields that 


u. = .NY 
i l- - 


(2.17) 


where .N ; i = 1, 2, 3 denote the three rows of [N]. 

l- 

In (2.17), since Y are not spatially dependent, we see that 


3u. 

l 


•3(5 1 ,C 2 ,C 3 ) } ~ [ 3(5 1 ,C 9 ,C,) v i 


(,N T )]Y 


3 2 u. 


’l ,s 2’*3 

.2 


{ 


r> - [- 


3(C 1 ,C 2 .C 3 ) 2/ l 3(C 1 ,5 2 ,5 3 ) 2 l “ 


(,n t )]y 


(2.18) 


(2.19) 


where here 


[ 3(5 lt 5 2 ,e 3 ) ( i- )] 


9 — ( N T ) 


- — ( n t ) 
3£ 2 i~ 

i — ( n t ) 

L 35 3 M- 


(2.20) 



7 


and 


[ d ~ 2 ( i NT)] 

9U r £ 2 ,5 3 ) 


r 





(.n t ) 


1- 


(.N T ) 


> 


V 


8 -(,n t ) 


9^3 1_ J 


( 2 . 21 ) 


The details of the differentiation of the actual components of [N] 

S a 2 3 2 

is given in Tables (2.2) - (2.4). Here j£—( )» 2 ( ^ an ^ ^ are 

depicted. Based on these expressions, the various other derivative ex- 
pressions can be generated. 

In the context of the foregoing nomenclature, the following expression 
can be generated for 3u^/a(n^,n 2 , n^) namely 


3 u i 0 

{ a(n lt n 2 ,n 3 j } " [J] 


r(.N T )]Y 


(2.22) 


where ie[l,3]. In a similar manner, (2.14), (2.20) and (2.21) yield the 
relation 

f — | s 

1 a(n 1 ,n 2 n 3 ) J 


[Jj] 


-1 


t — 2(.N T )]Y - [J3] 

3(s r c 2 ,s 3 ) 


[ 


3(C 1 »5 2 ,5 3 v i 


(,n t )]y 


(2.23) 


For transient situations, Part I derived the following FE formulation 


namely 
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/ { [G] T [S] [G] + [B*] T AS + [N] T ( [ mi ] [N] AY + [[m 2 HN] + 
R 

2[m 1 M[N])]Y + [tm 3 ][N] + [m 2 W[N]) + 


[m 1 ] Y ([N])]AY)}dv = AF 


(2.24) 


Upon use of the Newmark Beta type integration method, (2.24) was seen 
to reduce to the following more tractable form, that is 


t+at'tV + 'V + [ V ] “ ■ t + Pt ( S U a)S + a - + V 

where here 


(2.25) 


t+At 


V 


/ [N] T [C-[m 1 ] [N] + c ([m,][N] + 

R 31 ° 2 

2[m 2 ]T([N])) + [m 3 )[N] + [n^MlN]) + [n^M [N] )]dv 


(2.26) 


As can be seen again, consistent mass terms are embedded in (2.25) which 
involve the 'P( ) and y( ) differential operators. These can be treated 
through the use of (2.22) and (2.23) wherein the various derivatives of 
the components of [N] can be obtained/generated from Tables 2. 1-2.4. 


3. Thick Shell Element Formulation 

Certain engineered structures have the property that their normal 
strain components in the thickness direction are negligible in comparison 
to the in-plane values. In these cases, the structure can be modelled as 
a shell. For such situations, enhanced computational savings are obtained 
since the normal components in the thickness direction are entirely neglected. 
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I 


To allow for this modelling feature in steady and transient moving load 
problems, an eight node isoparametric shell element is developed. It is 
obtained by degenerating a three dimensional solid element [ 5 , 6 ]. A 
typical shell element and its associated nodes is shown in Fig. 3 . 1 . 

To start the development, it is noted that the in-plane shape function 
components are the same as those for the pl ane (£3=0) in t * ie s °l id 

element. The relative in-plane displacements between the top and bottom 
surfaces of the solid element are handled by rotational degrees of freedom 
and the thickness. These rotational degrees of freedom are with respect to 
the local axes, as shown in Fig. 3 . 2 . The element has three displacement 
and two rotational degrees of freedom at each node. 

The displacement field in the element, for the i th global direction is 


expressed as 


8 


u i - ° l N k <5 l ,C 2 )u i + ? N k (t l’ { 2 )h k ( ' V 2i 6 l +, 'h 6 2 ) 


( 3 . 1 ) 


k=l 


where u^ is the displacement of node k in the i direction, h^ is the 
thickness at the given node, (v!^, v!^) are components of the local vector 
in the (1,2) directions and lastly (6*, 6*) are the rotational degrees of 
freedom about the local ( 1 , 2 ) axes at the k th node. The derivatives with 
respect to the in-plane variables ^ and $ 2 of ( 3 * 1 > wil1 be of the same 
form, except that the shape function is replaced by the appropriate deri- 
vative. For example, the displacement derivative with respect to £3 is 


du. 
1 


8 

I ( 

k=l 


3N k (C l* C 2 ) 


k 


^ir^i k k k kw 

+ r at? - h k<-2i e ? + v ii e 2 )) 


( 3 . 2 ) 


! 
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Similar expressions can be derived for the remaining first and second 
derivatives. Note, the values of the requisite components of the shape 
function and their associated derivatives can be found using Tables (2.1) 
( 2 . 4 ) by setting £3 to zero. 

The first derivative of the displacement components with respect to 
the out of plane coordinate £^ is 8 iven b Y 



Z 2 N k^l* C 2 )h k^‘ V 2 i 6 l + v li 6 2 ) 
k=l 


(3.3) 


2 2 

The second derivative namely 3 u^/3£^ is zero. Continuing, the second 
derivatives with respect to (£^,£ 3 ) and (£ 2 *^ 3 ) ere a combination of (3.2) 
and (3.3). The forms of (3.2) and (3.3) and the appropriate values of the 
shape functions and their derivatives are used to give the quantities needed 
to define the transformed inertia matrices. 

The calculation of the transformed mass matrix using the formulation 
just described requires extensive numerical integration. Obviously, it 
would be preferred from a computational point of view to have a closed form 
expression for the transformed mass matrix. This is possible for the special 
case of rotation about a single axis. 

As a simplification of the formulation, considering body of revolution 
coordinates, the inertia in the circumferential and radial directions can 
be handled by the consistent approach [7], while the meridional direction 
inertia is handled by the lumped approach [7]. This is chosen because the 
inertia due to rotation lies in the circumferential - radial plane. Also, 
this way of expressing the problem will simplify the formulation and allow 
for closed form integration of the transformed mass matrix. Shape functions 
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for the consistent part of the formulation are written with respect to the 
elements local coordinate system rather than the isoparametric system. 

This puts a restriction on the shape and orientation of the element - its 
sides must lie in the circumferential, meridional, and radial directions. 

This is really not a restriction in body of revolution analyses since such 
an approach is the easiest way to build the finite element grid on a specific 
geometry. 

Overall, the foregoing simplifications can be implemented in several 
steps namely: 

i) Circumferential direction treated consistently; 

ii) Radial and meridional directions treated in lumped 
parameter format. 

Based on such an approach, the calculation of the matrix involves a 
quadratic variation in displacement in the circumferential direction 
specifically 

u*a + bn + cn^ (3.4) 

such that n is the circumferential coordinate. Note these proportionalities 
are assumed for the inertia calculations only. The displacement fields used 
for the stiffness matrix calculation remain as orginally defined. 

Using the known circumferential positions of the nodes, as shown for one 
circumferential line of nodes in Figure 3.2, the values of the coefficients 
in equation (3.4) are determined. These are then used to give shape functions 
and derivatives of the shape functions for these nodes. Each circumferential 
line of nodes in the element is defined as a zone, as shown in Figure 3.3. 
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The center zone has only end nodes in the shell element, but a pseudo-node 
is defined at the element center so that the same formulation can be used 
for this zone. Each zone is given a weighting value due to the lumped formu- 
lation being used in the meridional direction. 

The mass terms for the center node of the middle zone, the pseudo-node, 
must be condensed out because it is not present in the stiffness matrix 
calculations, and hence has no degrees of freedom in the global equations. 
This can be done in the usual manner on the element level, before the element 
mass matrix is assembled into the global mass matrix. This method of cal- 
culating the transformed mass matrix should be computationally more efficient 
because numerical integration is not needed. 

4. Three Dimensional Pantographing Gap Element Strategy 

There are several different ways to handle the problem of contact. 

These include: 

i) Direct application by continuous modification 
of boundary conditions [8]; 

ii) Hughes type auxiliary matrices [9]; 

iii) Penalty methods [10]} 

iv) Influence coefficients [11] and; 

v) Gap type procedures [12,15]. 

Each of these schemes have various advantages and disadvantages. Perhaps 
the single most often occurring problem involves the fact that some form 
of stiffness update and inversion is typically required during the overall 
solution process. 

To reduce computational effort, use is made here of substructuring 



concepts. Specifically, the global incremental stiffness formulation is 
recast in the form 
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(4.1) 


such that I, C, (AY^, AY C >, (AFj, AF c ) and ( [Kj] , . . . , [K c ] ) respectively 
define the subscripts associated with internal and contact degrees of 
freedom, the incremental deflection fields, the incremental nodal forces 
and lastly the various partitions of the tangential stiffness matrix. 
Solving for the incremental internal and contact degrees of freedom, we 
yield the following substructural expressions namely 

«C - [[K c ] - [K^uy-Vjt.naF,. - [KyiyVj (4.2) 

and 

AYj * -[K I ]" 1 [K IC ]AY C + [Kjl^AFj (4.3) 


To streamline the computational effort, [Kj] is updated only at the be - 
ginning of a given load step. During successive iterations, [K IC ], [ K C jl 
and [K r ] are intermittently updated depending on the type of contact pro- 

\j 

cedure employed. In this context, (4.2) and (4.3) can be recast in the 
following algorithmic form namely 


i A -C " [ i-l [i y ' i-l^CI 1 o tK I ] 1 i-l tK IC ]1 i A ?C 


-1 


i-l [K CI^ o [K I ] i-l*?I 


= o^I 3 " 1 i-l tK IC ] i A *C + o [K i rl i-l A ?I 


(4.4) 
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where the sub ( ) denotes the i th iteration of the given load step 
and sub-zero the initial value. 

To enable the more controlled handling of contact problems, local 
constraints can be imposed on the individual degrees of freedom associated 
with AY and AY-. Recalling Part 1 of the paper, the nonlinear FE formu- 

-I ~c 

lation can be solved via a multiply constrained partitioned NR solver. In 
terms of the substructuring noted in (4.1), the constraint process is applied 
directly to successive deflection excursions. Specifically, AY is replaced 
by [diag(A)jAY and hence (4.1) takes the form 


[Kj] 

w [K CI ] 


[K IC ] 

[K c ] 


[diag(A)] 




where 


[diag(A)] = 


~ [diag(Aj)] 

[ 0 ] 


[ 0 ] 

[diag(A c )] 


(4.5) 


(4.6) 


such that [diag(A I )] and [diag(A c )] respectively denote diagonal matrices 
defining individual constraints on the internal and contact degrees of 
freedom. Based on the form of (4.5), it follows that 


[diag( i A c )].AY c « 

WV - i-l [K CI ] i-l llt IC»i 4 !c - 

l-l [K Cl' o' K I> i-l 4 Jl 


and 
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tdiagC.Xj)]^ = 

o IK i rl i-l [K IC 1 [dia 8 < iV J i fi ?C + 


o'VVl 4 ?! 


(4.8) 


Note, each of the individual constraints appearing in [diag(^A^)] and 
[diag(^Aj)] can be defined by either: 

i) Establishing upper bounds on the deflections [3] of 
each particular substructure! partition of [Kj] and 
[K c ] or; 

ii) Employ a constraint function to bound allowable 
excursions [14]. 


For the problem of contact, it is preferable to employ i) wherein 
upperbound allowable excursions for the various partitions of AY^ can be 
continuously reset by monitoring the gap between contact surfaces. The 
various partitions of AY^ can be controlled in the manner defined by ii), 
namely by locally defined constraint functions [4]. 

Since a very large scale rolling simulation will be considered in this 
paper, to avoid the potential instabilities of the penalty method and in- 
efficiencies of the Hughes scheme, the pantographing gap methodology of 
Padovan and Moscarello [13] is employed. Noting Figs. 4.1, 4.2, and 4.3, 
pantographing enables the gap element to avoid severe distortion which 
occurs if ground nodes are constrained. Additionally, the gap element is 
subdivided into several zones associated with the level of integration em- 
ployed. Specifically, as can be seen from Fig. 4.4, since 3-3-3 integration 
is employed for the 20 node element, 9 separate zones are defined in the 



16 


contact area. These zones are alternatively turned on and off by the 
appropriate proximity criteria [13]. In much the way as elastic plastic 
formulations [ 15 ] . 

Noting Fig. 4.3, the element is continuously pantographed until a 
given node is found to satisfy the requisite proximity/contact criteria [13]. 
At this point, the zone associated with the given node is handled via an 
updated Lagrangian observer referencing the frozen pantographed portion 
of the node. If the node is found to release, then pantographing is resumed 
such that the prehistory is erased. 

The stiffness of the gap has essentially three phases of operation 
namely: 

i) Completely noncontacted wherein the stiffness is very 
small; 

ii) Partially or fully contacted wherein no slip is 
allowed and; 

iii) Partially or fully contacted wherein frictional slips 
is allowed. 

In the uncontacted mode, the gap element stiffness returned for assembly 
in the global matrix is numerical very small. For the partially or fully 
contacted mode wherein no slip is allowed, the appropriate zonal or whole 
element stiffness is made stiffer than the neighboring element. Lastly, for 
the case where frictional slip is allowed, the material stiffness of the gap 
is replaced by an adjusted orthotropic version of that of the neighboring 
structural element. Specifically: 

i) The direction normal to the contact surface is made very 
stiff (an order of magnitude larger than the neighboring 
element); 
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ii) The Poisson effect terms linking the normal and 
tangential behavior are deleted; 

iii) Tangential direction to contact is reset so as to have 
the same stiffness properties as neighboring element and; 

iv) Tangential contact nodal deflections are released and 
replaced by appropriate friction forces. 

Note, the criteria defining contact and slip-stick behavior are given 
by the following expressions namely: 


i) Contact conditions; 



- Y 


e 

£Cr it 


)-6 


f >0, no contact 
| <0 contact 


( 4 . 14 ) 


C >0 no contact 
\ <0 contact 


( 4 . 15 ) 


ii) Slip/stick conditions; 



F Crit> 


r >0 slip 
| <0 stick 


( 4 . 16 ) 


th © 6 

where here 2. defines the contact nodes of the e element, and Y^.^, 

5 p e p e an( i f respectively represent the normal deflection, contact 

* N2’ xl Crit 

gap distance of given node, error tolerance, normal force, tangential force 
and the slip force threshold. 

From an operational point of view. 


the slip inode will require the gap 
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stiffness to be restructured. In particular, given the e gap element 
expression, after partitioning we have that 



(4.9) 


where noting Fig. 4.4, the sub x denotes slipping tangential degrees of 
freedom at the contact nodes, while g defines the remaining field variables. 
Based on (4.9), the slip condition can be defined by the following sub- 
structural representation namely 

F® = [[K®] - [K® t ] [ K® ] ” 1 [ k®^ ] ] Y® + [K® t ][K®] _1 F t ' (4.10) 

Y® = -[K®] _1 [K® ]Y + tK®] -1 F e (4.11) 

-x x J xg -g x -x 

Depending on the number of element zones involved in frictional slipping, 
then 

F® = f®(F®) (4.12) 

- - 

where here denotes the normal forces at contact nodes and f the appro- 
-N ~ 

priately zonalized function defining the slip friction. In the classic 

0 Q 

case of Coulomb friction, the generalized function f (Fjj) is replaced 
by 

F® « [p®] F® (4.13) 

such that [y®] is properly structured so as to define the various zones 
of the element which are slipping. 
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As a final note, it follows that due to the use of moving coordinates, 
no special change needs to be implemented in the pantographing gap strategy. 
In particular it applies to both standing and rolling situations involving 
both steady and transient simulations. 

5. Benchmarking 

To benchmark the foregoing rolling 3-D shell and contact schemes, 
several comprehensive simulations will be considered. Overall, these 
include: 

i) 2-D evaluation to determine operational 
characteristics ; 

ii) Purely shell type 3-D model to evaluate mesh 
spacing requirements; 

iii) Mixed shell and 3-D element simulation 
of rolling tire. 

To enable such testing, the various formulations developed in sections 
2-4 were encoded in the NONSAP derivative code NFAP. The only major 
modification needed to convert the code involved the extension of the 
blocked skylined out of core solver [17] to handle the nonsymmetric 
matrices arising from the rolling/moving formulation. Beyond this, only 
standard modifications were introduced into the shell and 3-D (20 node) 
element libraries. 

Note, the main thrust of the benchmarking will be to: 

i) Correlate current modelling scheme with previous results 
and; 

ii) Compare results with experimental tests. 
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These comparisons will consider several aspects of system behavior, namely: 

i) Determine capabilities of rolling shell and 3-D 
elements to capture proper small deformation 
superposed on large eigenvalue properties and 
associated mode shapes; 

ii) Establish capabilities of gap contact scheme and; 

iii) Establish capabilities to handle overall interplay 
between geometric/material nonlinearity and inertia 
effects. 

As a first test, we consider the dynamic response of the 2-D ring on 
elastic foundation model of the tire [17]. Noting the model defined in 
Fig. 5.1, Figs. 5.2 and 5.3 illustrate the mixed shell/3-D and purely 3-D 
simulations. To evaluate and compare the dynamic characteristics of the 
shell and 3-D elements, we shall consider the problem of a circumferential 
traveling radial load. Fig. 5-4. Recalling the analytical solution developed 
by Padovan [18], when the circumferential traveling speed of the radial load 
is matched with the ratio of an individual frequency and its mode number 
M, namely uj^/M, a resonance type response is excited. For instance, con- 
sidering the 11 th mode depicted in Fig. 5.5, the shell and 3-D models yielded 
critical speeds of 119.7 mph and 120.3 mph respectively. Similar levels of 
accuracy were noted for all the critical speeds /natural frequencies ranging 
from the minimum to those above 180 mph (namely the first 25 frequencies). 

As a next test, we shall consider the problem of viscoelastic rolling contact 
Figure 5.6 illustrates the gap element supported model. Based on this simu- 
lation, Fig. 5.7 depicts the rolling contact shape and associated normal pres 
sure distribution in the contact zone. Both the shell and 3-D simulations 
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yielded essentially the same results for all reasonable ranges of rolling 
speed . 

The next model considered consists of simulating a two layer pressu- 
rized torus whose overall geometry and material properties are depicted 
in Fig. 5.8. As with the ring model, the pressurized torus is subject to a 
circumferentially traveling radial load. Figure 5.9 illustrates the changes 
in crown mode shape as a critical speed is approached and passed. The overall 
resonance mode shape is given in Fig. 5.10. These results were correlated 
with frequencies defined by the classic eigenvalue type formulation of the 
problem using a highly refined model. This comparison enabled the cali- 
bration of the appropriate element sizing. For instance, it was found that 
for the given shell element type; 

i) Element arcs that subtended greater than 1/36 of the 
circumference tended to yield slow convergence and; 

ii) In the meriodional/cross sectional orientation, at 
least 18 elements were needed to yield adequate model 
resolution. 

With such element spacing s, all frequencies in the range of engineering 
interest showed at most 3% deviations between the traveling and classic 
eigenvalue formulations. Such sensitivity studies enabled the reduction of 
the overall size of the tire model which even so is quite extensive. 

As the culminating benchmark. Fig. 5.11 illustrates the 15,000 degree 
of freedom model of a tire. Overall, the model treated the various internal 
laminations via Halpin Tsai [19] type correlations. The reason for the level 
of both meridional and circumferential refinement follows from the torus tests. 
These revealed that fairly uniform circumferential element spacing is needed 
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to define the proper mode shape/dynamic characteristics. 

The first test applied to the model consisted of pressurization and 
subsequent loading into ground contact. Figure 5.12 illustrates the gap 
element model used to define the contact region. Based on this model, the 
pressurization and subsequent loading into ground contact defined the 
axle load deflection curve denoted in Fig. 5.13. As can be seen, the model 
correlated quite well with experimental data obtained from the Firestone 
Tire and Rubber Company. The actual deflected shape of the contact region 
of the model is shown in Fig. 5.14. 

Next, we shall consider the models' capability to capture the frequency 
characteristics of the tire. Again, noting Fig. 5.15, this is achieved by 
placing a circumferentially travelling radial load on the crown section. 
Figures 5.15-5.21 illustrate various aspects of the dynamics. Note, as we 
are strictly interested in defining inertia capturing capabilities, damping 
was deleted from this series of tests. 

For velocities lower than 90 mph, little inertial effects were excited 
by pure rolling. As an example of this. Fig. 5.15 illustrates the tire re- 
sponse at 90 mph. Under modest increases above 90 mph, resonances were 
noted. For instance. Fig. 5.16 illustrates the crown node response to the 
approach and passing of a resonance speed. As can be seen, the range of speed 
is quite tight. The amplitude behavior during the resonance passing process 
is denoted in Fig. 5.17. Note, without damping, the resonance response is 
unbounded. At the critical speed (peak amplitude) associated with the given 
mode. Figs. 5.18 and 5.19 illustrate the crown and overall tire behavior. 

These results were within 17 . error of the experimentally generated data. With 
further increases in speed, other modes were excited. Figures 5.20 and 5.21 
illustrate the crown and whole tire response behavior at 138 mph. 
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It is interesting to note that the first critical speed of a rolling structure 
prototypically involves higher order modes. As the speed is gradually raised 
from the first critical, both lower and higher order modes may be excited. 

This is clearly seen by comparing Figs. 5.19 and 5.21 which illustrate the 8 tl1 
(117 mph) and 6 th (138 mph) modal responses. 

Such behavior is a direct outgrowth of the fact that the critical modes 
are directly related to w^/M. Noting Fig. 5.22, we see that usually oj^ is 
monotone increasing in M. Based on this, uj^/M prototypically is nonmonotone 
for small M. Figure 5.23 illustrates such behavior. The dip in the curve 
illustrated, points to the fact that higher order modes can give rise to the 
first critical. As the speed is increased, noting the trends depicted in Fig. 
5.23, various lower and higher resonces are subsequently excited. Additionally, 
depending on the inherent curvature of the versus M behavior of the struc- 
ture, the critical velocities can either be tightly or loosely packed. Such 
trends are depicted in Fig. 5.24. 

When damping is introduced, the individual modes tend to merge into a 
single so-called standing wave which appears behind the contact region and 
attenuates in the circumferential direction [3]. This behavior is clearly 
seen in Fig. 5.7. As the speed is gradually increased, the wave length tends 
to elongate. Such behavior is an outgrowth of the dissipation of higher order 
modes • 

As the last benchmark problem, we shall consider the steady freely rolling 
tire. Since the upper bound behavior of rolling tires is the so-called stand 
ing wave problem, main emphasis will be given to model this form of dynamic 
response. Figure (5.25) illustrates the standing wave response developed on 
a road wheel type test rig. As can be seen, extremely large deformation cha- 
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racteristics are excited. Employing the model depicted in Fig. 5.11, the 
rolling behavior is obtained in several steps namely: 

i) Pressurize tire structure accounting for follower 
type forces; 

ii) Push rolling tire into contact by varying hub 
deflection incrementally or alternatively; 

iii) Push non-rolling tire into contact by varying 
the hub deflection; once converged, rolling 
velocity is incrementally increased. 

To evaluate the importance of friction to the global effects, both 
pure slip and stick conditions were considered. While this induced different 
very localized shear distributions in the tread area of the contact patch, 
no changes were recorded in the global dynamic resonse. In this context, a 
stick condition was employed for the freely rolling full scale model. This 
significantly reduces run times. 

Based on the foregoing. Figs. 5.26-5.29 illustrate the response under 
different ranges of viscoelasticity. As can be seen from Figs. 5.26-5.27, for 
the very lightly damped case, standing wave patterns appear both fore and 
aft of the contact patch. This follows from the fact that there was insuf- 
ficient damping to attenuate the inertial interplay in circumferential direction. 
Once sufficient damping is introduced, essentially all the fore interactions are 
attenuated. 

From a comparison of Figs. 5.25 and 5.29, we see that excellent accuracy 
is obtained. In this context, it is noted that the moving element approach 
if used inconjunction with the appropriate FE mesh spacing can handle the com- 
plex dynamics associated with real world traveling load problems. This includes 
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the possibility of handling: 

i) 2-D, 3-D and shell type formulations; 

ii) Contact with and without friction; 

iii) Small and large deformations effects and; 

iv) Viscoelastic behavior. 

6. Summary 

Based on the moving strategy developed in Part 1, this paper derived 
3-D, shell and contact algorithm extensions. These modelling capabilities 
were extensively benchmarked to evaluate their operational capabilities. 

This included both analytical and experimental correlations. As was seen, 
excellent agreement was obtained over a wide range of physical situations 
namely: 

i) Traveling load problems involving moving 
velocities over the full interval including 
critical speeds; 

ii) 3-D, shell and 2-D situations; and 

iii) Full simulation of rolling contact accomo- 
dating viscoelasticity and the full definition 
of inertial effects. 

Note, due to the manner of formulation, the moving element procedure can 
be encoded into any of the currently available general purpose codes. This 
will enable the handling of moving load/boundary condition problems. 
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Figure and Table Captions 


Table No. 


Caption 


2.1 


2.2 


2.3 


2.4 


Shape function components for 20 node isoparametric 
solid element 

First derivative of shape function components with 
respect to 20 node isoparametric solid element 
Second derivative of shape function components with 
respect to 20 node isoparametric solid element 
Mixed second derivative of shape function components 
with respect to $ 2 ; 20 node isoparametric solid 


element 
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Figure and Table Captions 


Fig. No. 


2.1 


2.2 


3.1 


3.2 


3.3 


4.1 


4.2 


4.3 


4.4 


5.1 


5.2 


5.3 


5.4 


5.5 


5.6 


5.7 


5.8 


5.9 


Caption 

Three Dimensional 20 Node Solid Element in Cartesian Space 

Three Dimensional 20 Node Solid Element in £i» £2* £3 
Isoparametric Space 

Reference Surface of the Thick Shell Element in Cartesian 
Space 

Local Circumferential Positions for a Line of Nodes Used 
in the Closed Form Mass Formulation 

Division of the Shell Element into Three Zones for the 
Closed Form Approach 

Attachment and Constrained Face of the Three Dimensional 
Gap Element with Node Ordering 

Deformed Shape of a Gap Element with Lower Nodes Constrained 

Initial (- * -) and Pantographed ( — * — ) Shapes Minimizing 
Gap Element Distortion 

Zone Divisions in the Gap Element Used to Impose Partial 
Contact by Material Property Change at Integration Points 

Ring on Elastic Foundation Tire Model 

Finite Element Grid of the Ring on Elastic Foundation Tire 
Model Using Shell Elements for the Ring 

Finite Element Grid of the Ring on Elastic Foundation Tire 
Model Using Solid Elements for the Ring 

Line Load Across the Width of the Model Used as "Point" 

Load Excitation 

11th Mode Resonance Response at 145.6 Rad/Sec (119.7 MPH) 
for Ring Modeled by Shell Elements 

Finite Element Grid for the Ring on Elastic Foundation Model 
Including Gap Elements 

Deformed Shape and Contact Stress Distribution at 145.6 
Rad/Sec (119.7 MPH) and Damping (y=10 - *) 

Finite Element Grid of the Half Torus Model Viewed at a 
Slight Angle from the Axis of Rotation 

Response of the Torus Model as a Resonance/Critical Speed 
is Approached and Passed 



Fig. No. 

5.10 

5.11 

5.12 

5.13 

5.14 

5.15 

5.16 

5.17 

5.18 

5.19 

5.20 

5.21 

5.22 

5.23 

5.24 

5.25 


Caption 

4th Mode Resonance Response of the Torus Model at 284 
Rad/Sec 

Finite Element Grid of the Half Tire Model Used in the 
Rolling Contact Model (15000 Degrees) 

Finite Element Grid of the Quarter Tire Model with Gap 
Elements Attached in Potential Contact Zone 

Predicted and Measured Static Load-Deflection Response of 
Tire 

Deformed Shape of the Quarter Tire Model with 1.0 Inch 
Axle Deflection 

Side View of the FE Model Response to Radial Circumfer- 
entially Traveling Load Moving at 135.2 Rad/Sec (90 MPH) 

Predicted Response of Tire's Crown Nodes as Traveling 
Load Speed Approaches and Passes a Resonance/Critical 
Speed 

Traveling Speed-Radial Displacement Spectrum; Approach 
and Passing Behavior about Resonance/Critical Speed 

Response of the Tire's Crown Nodes at 175.7 Rad/Sec 
(117 MPH) Due to Circumferentially Traveling Radial Load 

Side View of the FE Model Response at 175.7 (117 MPH) 

Due to Circumferentially Traveling Radial Load 

Response of the Tire's Crown Nodes at 207 Rad/Sec (138 
MPH) Due to Circumferentially Traveling Radial Load 

Side View of the FE Model Response at 207 Rad/Sec (138 
MPH) Due to Circumferentially Traveling Radial Load 

Typical Frequency Spectrum Characteristics to Circum- 
ferential Mode Variations 

Critical Velocity Relationship 

Effects of Frequency Spectral Properties on Critical 
Speed Characteristics 

Standing Wave Response of Rolling Tire 

Response of the Damped Tire Model (y=10~^) Rotating at 
175.7 Rad/Sec (117 MPH) in Contact with the Ground Due 
to a .1 Inch Axle Deflection 


5.26 



Fig. No. 

5.27 

5.28 


Caption 

Response of Damped Tire Model (p“10 - ^) Rotating at 175.7 
Rad/Sec (117 MPH) in Contact with Ground Due to 1. Inch 
Axle Deflection 

Response of Damped Tire Model (y=1.5*10 - ^) Rotating at 
175.7 Rad/Sec (117 MPH) in Contact With the Ground Due 
to a 0.1 Inch Axle Deflection 

Response of the Damped Tire Model (y=1.5*10 - ^) Rotating 
at 175.7 Rad/Sec (117 MPH) in Contact with the Ground Due 
to a 1. Inch Axle Deflection 
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NODE 


(5 1 ; S 2 ; C 3 ) = (r; s; t) 

SHAPE FUNCTION (N) 


1 

(1/8) (1+r) (1+s) (1+t) - 

(1/2)(N 9 +N 12 +N 17 ) 

2 

( 1/8 ) { 1-r) ( 1+s) ( 1+t) - 

(1/2)(N 9 +N 10 +N 18 ) 

3 

(1/8) (1-r) (1-s) (1+t) - 

(i/2)(n 10 +n u +n 19 : 

4 

(1/8) (1+r) (1-s) (1+t) - 

(1/2)(N u +N 12 +N 20 

5 

(1/8) (1+r) (1+s) ( 1-t) - 

(1/2)(N u +N 16 +N 17 

6 

( 1/8) ( 1-r) ( 1+s) (1-t) - 

(1/2XN 13 +N 14 +N 18 

7 

(1/8) (1-r) (1-s) (1-t) - 

(1/2)(N 14 +N 15 +N 19 

8 

(1/8) (1+r) (1-s) (1-t) - 

(1/2)(N 15 +N 16 +N 20 

9 

(1/4) ( 1-r 2 ) (1+s) ( 1+t) 


10 

(1/4) (1-r) (1-s 2 ) ( 1+t) 


11 

(1/4) (1-r 2 ) (1-s) (1+t) 


12 

(1/4) (1+r) (1-s 2 ) ( 1+t) 


13 

( 1/4) ( 1-r 2 ) (1+s) ( 1-t) 


14 

(1/4) (1-r) (1-s 2 ) (1-t) 


15 

(1/4) ( 1-r 2 ) (1-s) ( 1-t) 


16 

(1/4) (1+r) (1-s 2 ) (1-t) 


17 

(1/4) (1+r) (1+s) (1-t 2 ) 


18 

(1/4) (1-r) (1+s) (1-t 2 ) 


19 

(1/4) (1-r) (1-s) (1-t 2 ) 


20 

(1/4) (1+r) (1-s) (1-t 2 ) 


Table 2.1 

Shape function components for 20 node 
isoparametric solid element 



NODE 

DERIVATIVE OF SHAPE FUNCTION (N»_) 

1 

( 1/8 ) ( 1+s) ( 1+t) - (1/2)(N 9 #r +N l2 #r *N 1? r ) 

2 

(-1/8) (1+s) (1+t) - (l/2)(N 9#r +N 10 #r + N 18 r ) 

3 

(-1/8) ( 1-s) (1+t) - (l/ 2 HN 10#r +N 1Ur +N 19 r ) 

' 4 

(1/8) (1-s) (1+t) - (l/2)(N n /r + N 12#r +N 20 r ) 

5 

(1/8) (1+s) (1-t) - (l/ 2 )(N 13#r +N 16#r +N 17#r ) 

6 

(-1/8) (1+s) (1-t) - (1/2)(N 13 r +N 14 fr +N 18 v ) 

7 

(-1/8) (1-s) (1-t) - < 1 / 2 )(N 14 r +N 15 r +N 19 r ) 

8 

(l/8)(l-s) (1-t) - (1/2)(N 15 #r +N lg #r +N 20#r ) 

9 

(-1/2) (r) ( 1+s) (1+t) 

10 

(-l/4)(l-s 2 )(l+t) 

11 

(-1/2) (r) ( 1-s) (1+t) 

12 

( 1/4) ( 1-s 2 ) (1+t) 

13 

(-1/2) (r) (1+s) (1-t) 

14 

(-1/4) ( 1-s 2 ) ( 1-t) 

15 

(-1/2) (r) (1-s) (1-t) 

16 

( 1/4) ( 1-s 2 ) (1-t) 

17 

( 1/4) { 1+s) ( 1-t 2 ) 

18 

(-1/4) ( 1+s) (1-t 2 ) 

19 

(-1/4) (1-s) (1-t 2 ) 

20 

(1/4) (1-s) (1-t 2 ) 

Table 2.2 

First derivative of shape function components 


with respect to 20 node isoparametric solid 
element 


NODE 


DERIVATIVE OF SHAPE FUNCTION (_Nf_ ri J 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


<-l/2HN 9(rr M! 12(rr *N 17iI . r > 

(-1/2) ( N 9 #r c +N 10,rr +N 18 ,rr* 

(-1/2) (N 10 , rr +li ii ( rr + **19.rr^ 

1 - 1/2 ) < N 1 1 , rr +H i 2 , tr +N 20 , r r* 

(-1/2) (**i 3 ( rr* N 16#rr + **17 ,cr* 

(-l/2)(N l j rr +H l4>rr »N l8>rt ) 

(-l/2)(N 14jCt »N l5>rr +N 19<rr ) 

(-l/2)(Ni 5irr ^i 6< rr tN 20.rr' 

(-1/2) ( 1+s) (1+t) 

0 

(-1/2) ( 1-s) (1+t) 

0 

(-1/2) (l+s)(l-t) 

0 

(-1/2) ( 1-s) ( 1-t) 

0 

0 

0 

0 

0 

Second derivative of shape function components 
with respect to ; 20 node isoparametric solid 
element 


Table 2.3 



NODE 


(5-^ £ 2 * ^3^ E ( r; s; 

DERIVATIVE OF SHAPE FUNCTION (N, rj ) 


1 

(1/8) (1+t) - (l/2HN, trs + »i2,rs +N 17,rs > 

2 

(-1/8) (1+t) - (l/JHN 9>rs +N 10>rs *N 18>rs ) 

3 

(1/8) (1+t) - U/2)<N 10 ,rs +,, ll.ts +N 19,rs ) 

4 

(“1/8) (1+t) - (1/2) ^ fr 3+**l2 # rs + **20,rs^ 

5 

(1/8) (1-t) - (l/2)(N l3trJ +H 1<<rs +l' 17 , rs » 

6 

(“1/8) (l“t) - (l/2)(N 1 j <rl +H l4>tl +>l 18 , rs > 

7 

( 1/8) ( 1— t) - (l/2)(N| 4 ^ rJ +N 3 j^ rs +N 39 ^ r# ) 

8 

(-l/8)(l-t) - (l/2)(N l5iCS +N Uirs +N 20iI . s ) 

9 

(-1/2) (r) (1+t) 

10 

( 1/2) (s) ( 1+t) 

11 

( 1/2) ( r) ( 1+t) 

12 

(-1/2) (s) (1+t) 

13 

(-1/2) (r) (1-t) 

14 

( 1/2) (s) ( 1-t) 

15 

(1/2) (r) (1-t) 

16 

(-1/2) (s) (1-t) 

17 

(1/4) (1-t 2 ) 

18 

(-1/4) (1-t 2 ) 

19 

(1/4) (1-t 2 ) 

20 

(-1/4) (1-t 2 ) 

Table 2.4 

Mixed second derivative of shape function components 
with respect to 20 node isoparametric solid 

element 
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Fig . 3.3 

division of the shell 

ELEMENT INTO THREE ZONES 
FOR THE CLOSED FORM 
APPROACH 



12 
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16 
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18 


19 


CONSTRAINED 


A . 1 ATTACHMENT AND CONSTRAINED 
FACES OF THE THREE DIMEN- 
SIONAL GAR ELEMENT WITH NODE 
ORDERING 


z 


A 



CONSTRAINED 


Fig- A. 2 DEFORMED SHARE OF A GAR 

ELEMENT WITH LOWER NODES 
CONSTRAINED 



INITIAL C — — > AND 

PANTOGRAPHED C 3C - ) 
SHAPES MINIMIZING GAP 
ELEMENT DISTORTION 



ZONE DIVISIONS IN THE GAP 
ELEMENT USED TO IMPOSE 
PARTIAL CONTACT BY MATERIAL 
PROPERTY CHANGE AT INTE- 
GRATION POINTS 




3-D SOLID 
ELEMENT 
< FOUNDATION ) 




SHELL ELEMENT 
(RING) 


Fig . 5.2 FINITE ELEMENT GRID OF THE 

RING ON ELASTIC FOUNDATION 
TIRE MODEL USING SHELL 

ELEMENTS FOR THE RING 




3-D SOLID 
ELEMENTS 
C FOUNDATION ) 


Fig . 


\ 

3-D SOLID ELEMENTS 
(RING) 


5.3 FINITE ELEMENT GRID OF THE RING 

ON ELASTIC FOUNDATION TIRE MODEL 
USING SOLID ELEMENTS FOR THE 
RING 




APPLIED "TRAVELING" 
RADIAL "POINT" LOAD 


Rig. 5. A LINE LOAD ACROSS THE WIDTH 

OF THE MODEL USED AS "POINT" 
LOAD EXCITATION 





5.5 llth MODE RESONANCE 

RESPONSE AT 145 . €> RAD/ SEC 
(119.7 MPH) FOR RING 
MODELED BY SHELL ELEMENTS 


GAP ELEMENTS 



FINITE ELEMENT GRID FOR 
THE RING ON ELASTIC 
FOUNDATION MODEL INCLUDING 
GAP ELEMENTS 




Fig . 5.7 DEFORMED SHAPE AND CONTACT STRESS 

DISTRIBUTION AT 145 . 6 RAD/ SEC 
(119.7 MPH) AND DAMPING 
(y=10“ 4 ) 



I 


RADIAL TRAVELING LOAD 

Fig . 5.8 FINITE ELEMENT GRID OF THE HALF 

TORUS MODEL VIEWED AT A SLIGHT 
ANGLE FROM THE AXIS OF ROTATION 





F ±g . 5.11 FINITE ELEMENT GRID OF THE 

HALF TIRE MODEL USED IN THE 
ROLLING CONTACT MODEL 
(15000 DEGREES) 



GAP ELEMENTS 

Fig . 5.12 FINITE ELEMENT GRID OF THE 

QUARTER TIRE MODEL WITH GAR 
ELEMENTS ATTACHED IN 

POTENTIAL CONTACT ZONE 


LOAD ( LB ) 



DEFLECTION ( INCHES > 


Fig. 5.13 PREDICTED AND MEASURED 

STATIC LOAD -DEFLECTION 
RESPONSE OF TIRE 


Fig. 5.14 DEFORMED SHAPE OF THE 

QUARTER TIRE MODEL WITH 
1 . O INCH AXLE DEFLECTION 



| RADIAL TRAVELING 
LOAD 

F± s . 5.15 SIDE VIEW OF THE FE 

MODEL RESPONSE TO RADIAL 
CIRCUMFERENTIALLY TRAVEL- 
ING LOAD MOVING AT 135 . 2 
RAD/SEC (90 MPK) 


O RAD/SEC f ) 176.5 RAD/SEi 

pred icted response op tire 

CROWN NODES AS TRAVELING L 
SPEED APPROACHES AND PASSE 
RESONANCE/ CRITICAL SPEED 




170 180 
ROTATIONAL. SPEED (RAD/SEC) 


F±g. 5.17 TRAVELING SPEED-RADIAL 

DISPLACEMENT SPECTRUM i 
APPROACH AND PASSING 
BEHAVIOR ABOUT RESONANCE/ 
CRITICAL SPEED 


| RADIAL TRAVELING 
LOAD 

RESPONSE OF THE TIRE'S 
CROWN NODES AT 175.7 
RAD/ SEC (117 MPH) DOE TO 
CIRCUMFERENTIALLY TRAVELING 
RADIAL LOAD 




LOAD 

Fig. 5.19 SIDE VIEW OF THE FE MODEL 

RESPONSE AT 175.7 (117 MPH) 

DUE TO CIRCUMFERENTIALLY 
TRAVELING RADIAL LOAD 



• RADIAL TRAVELING 
LOAD 

Fig. 5.20 RESPONSE OF THE TIRE * S CROWN 

NODES AT 207 RAD/SEC (138 MPH ) 
DUE TO CIRCUMFERENTIALLY 
TRAVELING RADIAL LOAD 



Fig . 



LOAD 

5.21 SIDE VIEW OF THE FE MODEL 
RESPONSE AT 207 RAD/SEC 
<138 MPH) DUE TO CIRCUMFER- 
ENTIALLY TRAVELING RADIAL 
LOAD 
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Fig 
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5.22 TYPICAL FREQUENCY S PECTRUM 
CHARACTERISTICS TO CIRCUM- 
FERENTIAL MODE VARIATIONS 



CRITICAL ^ co m /M ^ 







Fig- 5.26 RES PONSE OF THE DAMPED TXRE 

MODEL (y=10" 4 ) ROTATING AT 
175.7 RAD/SEC (117 MPH) IN 
CONTACT WITH THE GROUND DUE 
TO A . 1 INCH AXLE DEFLECTION 



±g. 5.27 RESPONSE OP DAMPED TIRE MODEL 

C vi=10 “ A ) ROTATING AT 175-7 
RAD/SEC (117 MPH) IN CONTACT 
WITH GROUND DUE TO 1 . INCH 
AXLE DEFLECTION 



Fig . 5-28 RESPONSE OF DAMPED TIRE MODEL 

C u= 1.5*10 _4 ) ROTATING AT 175.7 
RAD/ SEC (117 MPH) IN CONTACT 
WITH THE GROUND DUE TO A. 0.1 
INCH AXLE DEFLECTION 


STANDING 

WAVE 


RESPONSE OF THE DAMPED TIRE 
MODEL ( vi= 1.5*10" 4 ) ROTATING 
AT 175.7 RAD/SEC (117 MPH) 

IN CONTACT WITH THE GROUND 
DUE TO A 1 . INCH AXLE DEFLEC 
TION 


